#define _REAL_ double precision

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!+ [Enter a one-line description of subroutine grtopr here]
subroutine grtopr(x,y,z,xx,yy,zz,i0,j0,k0,qg,inf1,inf2,inf3, &
      q0,qx,qy,qz,qxx,qyy,qzz,qxy,qxz,qyz)
   implicit _REAL_ (a-h,o-z)
   parameter (nq = 27)
   
   !       **************************************************************
   !       *                                                            *
   !       *  GrToPr find:                                              *
   !       *     if inf1=1: the value q0 of a quantity q at (xx,yy,zz)  *
   !       *     if inf2=1: the value qx,qy,qz of  q at (xx,yy,zz)      *
   !       *     if inf3=1: the value qxx,qyy,qzz,qxy,qxz,qyz of q at   *
   !       *                (xx,yy,zz)                                  *
   !       *  based upon its values at grid points qg(*,*,*).           *
   !       *                                                            *
   !       *   im=1: 8-pint bilinear interpolation scheme used.         *
   !       *         (only used in case: inf1=1 only)                   *
   !       *   im=2: General interpolation scheme used.                 *
   !       *                                                            *
   !       **************************************************************

   common /lmn/l, m, n, nirreg
   common /hxyz/hx,hy,hz,hmax

   dimension x(0:l+1),y(0:m+1),z(0:n+1)
   dimension qg(0:l+1,0:m+1,0:n+1)

   dimension w1(10, nq), w4(10)
   dimension sd(11), uw(10,10), v(nq, nq)
   dimension ew(nq)
   dimension ew1(nq), ew2(nq), ew3(nq)
   dimension ew4(nq), ew5(nq), ew6(nq)
   dimension ew7(nq), ew8(nq), ew9(nq), ew10(nq)
   dimension w31(nq), w32(nq), w33(nq)
   dimension w34(nq), w35(nq), w36(nq)
   dimension w37(nq), w38(nq), w39(nq), w310(nq)
   dimension qcoe(10,nq)

   dimension sd2(10), work(1000)

   ! ----- set SVD subroutine: isvd=1  call ssvdc()
   !                           isvd=2  call sgecvd()

   isvd = 1

   im = 2

   !       i0 = nint((xx - x(0))/hx)
   !       j0 = nint((yy - y(0))/hy)
   !       k0 = nint((zz - z(0))/hz)


   if (im == 1 .and. inf2 == 0 .and. inf3 == 0) then
      i1 = i0 + 1
      j1 = j0 + 1
      k1 = k0 + 1
      x0 = x(i0)
      y0 = y(j0)
      z0 = z(k0)
      x1 = x0 + hx
      y1 = y0 + hy
      z1 = z0 + hz

      x0b = 2.0*(x1 - xx)/hx
      x1b = 2.0*(xx - x0)/hx
      y0b = 2.0*(y1 - yy)/hy
      y1b = 2.0*(yy - y0)/hy
      z0b = 2.0*(z1 - zz)/hz
      z1b = 2.0*(zz - z0)/hz

      q0 =   x0b * y0b * z0b * qg(i0,j0,k0) &
            + x0b * y0b * z1b * qg(i0,j0,k1) &
            + x0b * y1b * z0b * qg(i0,j1,k0) &
            + x0b * y1b * z1b * qg(i0,j1,k1) &
            + x1b * y0b * z0b * qg(i1,j0,k0) &
            + x1b * y0b * z1b * qg(i1,j0,k1) &
            + x1b * y1b * z0b * qg(i1,j1,k0) &
            + x1b * y1b * z1b * qg(i1,j1,k1)

      q0 = q0/8.0

   else

      do i=1, 10
         do j=1, nq
            qcoe(i,j) = 0.0
         end do
      end do

      ! -------- initialize

      do jj=1, nq
         do ii=1, nq
            v(ii, jj) = 0.0      !# the V matrix in SVD: U*S*V
         end do
         do ii=1, 10
            w1(ii, jj) = 0.0     !# the coefficient matrix to be SVDed
         end do
      end do

      ! -------- Generate the matrix

      nsub = 0
      do i1 = i0-1, i0+1
         do j1 = j0-1, j0+1
            do k1 = k0-1, k0+1
               if (i1 < 1 .or. j1 < 1 .or. k1 < 1 &
                     .or. i1 > l .or. j1 > m .or. k1 > n) goto 55

               nsub = nsub + 1
               w1(1, nsub)  = 1.0
               w1(2, nsub)  = x(i1) - xx
               w1(3, nsub)  = y(j1) - yy
               w1(4, nsub)  = z(k1) - zz
               w1(5, nsub)  = 0.5*(x(i1) - xx)*(x(i1) - xx)
               w1(6, nsub)  = 0.5*(y(j1) - yy)*(y(j1) - yy)
               w1(7, nsub)  = 0.5*(z(k1) - zz)*(z(k1) - zz)
               w1(8, nsub)  = (x(i1) - xx)*(y(j1) - yy)
               w1(9, nsub)  = (x(i1) - xx)*(z(k1) - zz)
               w1(10, nsub) = (y(j1) - yy)*(z(k1) - zz)
               55 continue
            end do
         end do
      end do

      if (nsub < 12) then
         write(*,*) "   nsub is too small in GrToPr()!"
         !             print *,xx,yy,zz
         stop
      end if

      ! -------- Call least square routine

      if (isvd == 1) then
         job = 11
         call dsvdc(w1,10,10,nsub,sd,ew,uw,10,v,nq,w4,job,inf)
      else
         call dgesvd('A','A',10,nsub,w1,10,sd2,uw,10,v,nq, &
               work, 1000, inf)
         do ij = 1, 10
            sd(ij) = sd2(ij)
         end do
      end if

      if (inf /= 0) then
         write(*,*) inf, " - ssvdc() or sgesvd() failed in GrToPr!"
         stop
      end if


      do i1=1, 10
         if (abs(sd(i1)) > 1.0e-14) then
            if (inf1 == 1) then
               ew1(i1)=uw(1, i1)/sd(i1)
            end if
            if (inf2 == 1) then
               ew2(i1)=uw(2, i1)/sd(i1)
               ew3(i1)=uw(3, i1)/sd(i1)
               ew4(i1)=uw(4, i1)/sd(i1)
            end if
            if (inf3 == 1) then
               ew5(i1)=uw(5, i1)/sd(i1)
               ew6(i1)=uw(6, i1)/sd(i1)
               ew7(i1)=uw(7, i1)/sd(i1)
               ew8(i1)=uw(8, i1)/sd(i1)
               ew9(i1)=uw(9, i1)/sd(i1)
               ew10(i1)=uw(10, i1)/sd(i1)
            end if
         else
            if (inf1 == 1) then
               ew1(i1)= 0.0
            end if
            if (inf2 == 1) then
               ew2(i1)= 0.0
               ew3(i1)= 0.0
               ew4(i1)= 0.0
            end if
            if (inf3 == 1) then
               ew5(i1)= 0.0
               ew6(i1)= 0.0
               ew7(i1)= 0.0
               ew8(i1)= 0.0
               ew9(i1)= 0.0
               ew10(i1)= 0.0
            end if
         end if  ! (abs(sd(i1)) > 1.0e-14)
      end do  !  i1=1, 10


      ! -------- w3(i) is the solution

      do i1=1, nsub
         if (inf1 == 1) then
            w31(i1) = 0.0
         end if
         if (inf2 == 1)  then
            w32(i1) = 0.0
            w33(i1) = 0.0
            w34(i1) = 0.0
         end if
         if (inf3 == 1) then
            w35(i1) = 0.0
            w36(i1) = 0.0
            w37(i1) = 0.0
            w38(i1) = 0.0
            w39(i1) = 0.0
            w310(i1) = 0.0
         end if

         do j1=1,10
            if (isvd == 1) then
               if (inf1 == 1) then
                  w31(i1) = w31(i1) + v(i1,j1)*ew1(j1)
               end if
               if (inf2 == 1) then
                  w32(i1) = w32(i1) + v(i1,j1)*ew2(j1)
                  w33(i1) = w33(i1) + v(i1,j1)*ew3(j1)
                  w34(i1) = w34(i1) + v(i1,j1)*ew4(j1)
               end if
               if (inf3 == 1) then
                  w35(i1) = w35(i1) + v(i1,j1)*ew5(j1)
                  w36(i1) = w36(i1) + v(i1,j1)*ew6(j1)
                  w37(i1) = w37(i1) + v(i1,j1)*ew7(j1)
                  w38(i1) = w38(i1) + v(i1,j1)*ew8(j1)
                  w39(i1) = w39(i1) + v(i1,j1)*ew9(j1)
                  w310(i1) = w310(i1) + v(i1,j1)*ew10(j1)
               end if
            else
               if (inf1 == 1) then
                  w31(i1) = w31(i1) + v(j1,i1)*ew1(j1)
               end if
               if (inf2 == 1) then
                  w32(i1) = w32(i1) + v(j1,i1)*ew2(j1)
                  w33(i1) = w33(i1) + v(j1,i1)*ew3(j1)
                  w34(i1) = w34(i1) + v(j1,i1)*ew4(j1)
               end if
               if (inf3 == 1) then
                  w35(i1) = w35(i1) + v(j1,i1)*ew5(j1)
                  w36(i1) = w36(i1) + v(j1,i1)*ew6(j1)
                  w37(i1) = w37(i1) + v(j1,i1)*ew7(j1)
                  w38(i1) = w38(i1) + v(j1,i1)*ew8(j1)
                  w39(i1) = w39(i1) + v(j1,i1)*ew9(j1)
                  w310(i1) = w310(i1) + v(j1,i1)*ew10(j1)
               end if
            end if  ! (isvd == 1)
         end do  !  j1=1,10

         if (inf1 == 1) then
            qcoe(1,i1) = w31(i1)
         end if
         if (inf2 == 1) then
            qcoe(2,i1) = w32(i1)
            qcoe(3,i1) = w33(i1)
            qcoe(4,i1) = w34(i1)
         end if
         if (inf3 == 1) then
            qcoe(5,i1) = w35(i1)
            qcoe(6,i1) = w36(i1)
            qcoe(7,i1) = w37(i1)
            qcoe(8,i1) = w38(i1)
            qcoe(9,i1) = w39(i1)
            qcoe(10,i1) = w310(i1)
         end if

      end do  !  i1=1, nsub

      q0  = 0
      qx  = 0
      qy  = 0
      qz  = 0
      qxx = 0
      qyy = 0
      qzz = 0
      qxy = 0
      qxz = 0
      qyz = 0

      nsub = 0
      do i1 = i0-1, i0+1
         do j1 = j0-1, j0+1
            do k1 = k0-1, k0+1
               if (i1 < 1 .or. j1 < 1 .or. k1 < 1 &
                     .or. i1 > l .or. j1 > m .or. k1 > n) goto 155
               if (qg(i1,j1,k1) > 999.0 ) then
                  print *,'unexcepted grid',i1,j1,k1,qg(i1,j1,k1)
                  print *,i0,j0,k0
                  print *,xx,yy,zz
                  stop
               end if

               nsub = nsub + 1
               if (inf1 == 1) then
                  q0  =  q0  + qcoe(1,nsub)*qg(i1,j1,k1)
               end if
               if (inf2 == 1) then
                  qx  =  qx  + qcoe(2,nsub)*qg(i1,j1,k1)
                  qy  =  qy  + qcoe(3,nsub)*qg(i1,j1,k1)
                  qz  =  qz  + qcoe(4,nsub)*qg(i1,j1,k1)
               end if
               if (inf3 == 1) then
                  qxx =  qxx + qcoe(5,nsub)*qg(i1,j1,k1)
                  qyy =  qyy + qcoe(6,nsub)*qg(i1,j1,k1)
                  qzz =  qzz + qcoe(7,nsub)*qg(i1,j1,k1)
                  qxy =  qxy + qcoe(8,nsub)*qg(i1,j1,k1)
                  qxz =  qxz + qcoe(9,nsub)*qg(i1,j1,k1)
                  qyz =  qyz + qcoe(10,nsub)*qg(i1,j1,k1)
               end if
               155 continue
            end do
         end do
      end do

   end if

   return
end subroutine grtopr 

